# Install Package to reverse geocoding. 
!pip install geopy

# Load Packages
import numpy as np
import pandas as pd 

# Load Business Location Data
sm = pd.read_csv(".../POIN_MASTER_010319.csv")

# pre-selected targetted list from the business location data.
 # in the following analysis, our group only focus on those stores.
chain_store_list = ['GEOCVS','GEODD','GEODICK','GEOHDPT','GEOHGOOD','GEOJCP',
                    'GEOKOHL','GEOKRGR','GEOLOWE','GEOMACYS','GEOMASS','GEOODEPT',
                    'GEOROSS','GEOSAFEW','GEOSAV','GEOSEAR','GEOTJ','GEOTRGT',
                    'GEOWALG','GEOWMT']

# remove unrelevant business from the original data frame. 
sm = sm[sm['gitext'].isin(chain_store_list)]

sm.info()

sm_isna = pd.read_csv(".../temp_sm_isna.csv")
sm_isnotna = sm[sm['zip'].isna() == False]

sm_isna.info()

sm_new = pd.concat([sm_isna,sm_isnotna],axis=0)
sm_new.drop(columns=['Unnamed: 0'],inplace=True)
# sm_new.to_csv("reverse_geocodone.csv")

# sm_new = pd.read_csv(".../reverse_geocodone.csv")
# sm_new.drop(columns=['Unnamed: 0'],inplace=True)

f = sm_new[['gitext', 'poiname', 'addr1', 'city', 'state','zip','latitude','longitude']]
# Remove canadian zipcode 
f = f[~f.zip.str.match('[a-zA-Z]')]

# select out the rows that zip is like xxxxx-xxxx
a = f[f.zip.str.match('^\d{5}-')]
a['zip'] = a['zip'].apply(lambda x:x[:5])

# select out the rows that zip is like xxxxx:xxxxx
b = f[f.zip.str.match('^\d{5}:\d{5}')]
b['zip'] = b['zip'].apply(lambda x:x[:5])

# select out the rows that zip length is less than 5 and add 0 in front of it
c = f[f['zip'].apply(lambda x: len(x)<5)]
c['zip'] = c['zip'].apply(lambda x: '{0:0>5}'.format(x))

# select out the rows that zip is like xxxxx:...
d = f[f.zip.str.match('^\d{5};')]
d['zip'] = d['zip'].apply(lambda x:x[:5])

# zipcode that is already in a desired format 
e = f[f.zip.str.match('^\d{5}$')]

# concat them together to get the final business location dataset 
final = pd.concat([a,b,c,d,e],axis=0)
# final.to_csv("full_zip.csv")

final.info()

# load the zipcodes from that are generated from the reverse-geocoding process.
fl = open(".../zip.txt",'w')
for i in zipz:
    fl.write(i+'\n')
fl.close()

zips = open('.../zip.txt', 'r').readlines()
zips = [z.replace('\n', '') for z in zips]
# print(len(zips))

# load the pre-selected attributes.
filed = open('.../field of interest.txt', 'r').readlines()
filed = [h.replace('\n', '') for h in filed]
print(len(filed))
filed = ','.join(filed)

df = pd.DataFrame(columns=["zipcode","B02001_001E","B02001_002E","B02001_003E","B02001_004E","B02001_005E",
                           "B02001_006E","B02001_007E","B02001_008E","B01003_001E","B08301_002E","B08301_003E",
                           "B08301_004E","B08301_005E","B08301_006E","B08301_007E","B08301_008E","B08301_009E",
                           "B08301_010E","B08301_011E","B08301_012E","B08301_013E","B08301_014E","B08301_015E",
                           "B08301_016E","B08301_017E","B08301_018E","B08301_019E","B08301_020E","B08301_021E",
                           "B08014_002E","B08014_003E","B08014_004E","B08014_005E","B08014_006E","B08014_007E"])

for i in range(102): 
  a = zips[i*100:(1+i)*100]
  a = ','.join(a)
  url = f'https://api.census.gov/data/2018/acs/acs5?key={apiKey}&get={filed}&for=zip%20code%20tabulation%20area:{a}'
  response = requests.get(url)
  #load the response into a JSON, ignoring the first element which is just field labels
  formattedResponse = json.loads(response.text)[1:]
  #flip the order of the response from [population, zipcode] -> [zipcode, population]
  formattedResponse = [item[::-1] for item in formattedResponse]
  #store the response in a dataframe
  zippp = pd.DataFrame(columns=["zipcode","B02001_001E","B02001_002E","B02001_003E","B02001_004E","B02001_005E",
                                "B02001_006E","B02001_007E","B02001_008E","B01003_001E","B08301_002E","B08301_003E",
                                "B08301_004E","B08301_005E","B08301_006E","B08301_007E","B08301_008E","B08301_009E",
                                "B08301_010E","B08301_011E","B08301_012E","B08301_013E","B08301_014E","B08301_015E",
                                "B08301_016E","B08301_017E","B08301_018E","B08301_019E","B08301_020E","B08301_021E",
                                "B08014_002E","B08014_003E","B08014_004E","B08014_005E","B08014_006E","B08014_007E"],
                       data=formattedResponse)
  df = pd.concat([df,zippp],axis=0)

df = df[df['B02001_001E'].isna() == False]

ma = pd.read_csv(".../mosaic_byZip.csv")

# format the mosaic data's zipcode column into standard xxxxx 
ma['zip'] = ma['zip'].apply(lambda x: '{0:0>5}'.format(x))

cs_and_ma = df.merge(ma,how='left',left_on='zipcode',right_on='zip')
cs_and_ma.drop(columns=['zip'],inplace=True)

# Creating wide table for further multi-labeled machine learning 
final = final[final['zip'].isin(zip_final)]
final = final[['gitext','zip']]
final = final.astype(str)
final['existence'] = [True] * len(final)
final_wide = final.pivot_table(index='zip',columns='gitext',values='existence')
final_wide.fillna(False,inplace = True)

# Combine the business location data with the census + mosaic dataset 
final = final_wide.reset_index()
final = final.merge(cs_and_ma,how='left',left_on='zip',right_on='zipcode')
final.drop(columns=['zipcode'],inplace=True)

# final.to_csv("fulldata.csv")

raw = final

# we drop columns where either do not contain ACS information or Mosaic information. 
raw.dropna(axis=0,inplace=True) # there are 9855 zipcodes left. 

# create our multi-lable dependent variables 
label = raw.loc[:,'GEOCVS':'GEOWMT'] 
y = label*1
# y.shape

# create our independent variabels into X
X = raw.loc[:,'B02001_001E':'political_affiliation_p111975_n']
X.shape

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_stan = scaler.fit_transform(X)

# Train test split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X_stan, y, test_size=0.2, random_state=101)

from sklearn.metrics import accuracy_score,roc_curve, auc, roc_auc_score
from sklearn.model_selection import KFold, cross_val_score
from sklearn import metrics
from sklearn.metrics import classification_report

n_folds = 10
def get_CVacc(model):
    """
    Return the accuracy score
    """
    # Set KFold to shuffle data before the split
    kf = KFold(n_folds, shuffle=True, random_state=42)
    
    # Get accuracy score
    accuracy_score = cross_val_score(model, X_stan, y, scoring="accuracy", cv=kf)
    
    return accuracy_score.mean()

def get_acc(model):
    """
    Return accuracy score
    """
    model.fit(X_train, y_train)
    predictions = model.predict(X_test)
    acc = accuracy_score(y_test,predictions)
    print(classification_report(y_test, predictions))
    return acc

from sklearn.ensemble import RandomForestClassifier

rfc = RandomForestClassifier(n_estimators=500)
print("accuracy_rfc: ", get_acc(rfc))

import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Activation,Dropout,BatchNormalization
from tensorflow.keras.constraints import max_norm

# Model 1:
model = Sequential()

# input layer
model.add(Dense(246,  activation='sigmoid'))
model.add(Dropout(0.5))

# hidden layer
model.add(Dense(64, activation='sigmoid'))
model.add(Dropout(0.5))

# hidden layer
model.add(Dense(128, activation='sigmoid'))
model.add(Dropout(0.5))

# hidden layer
model.add(Dense(32, activation='sigmoid'))
model.add(Dropout(0.5))

# output layer
model.add(Dense(units=20,activation='sigmoid'))

# Compile model
model.compile(loss='binary_crossentropy', optimizer='adam',metrics=['accuracy'])

# Model 2: 
model2 = Sequential()
model2.add(Dense(246, input_dim= 246, activation='sigmoid'))
model2.add(Dense(128, activation='sigmoid'))
model2.add(Dense(64, activation='sigmoid'))
model2.add(Dropout(0.5))
model2.add(Dense(64, activation='sigmoid'))
model2.add(Dense(32, activation='sigmoid'))
model2.add(Dropout(0.5))
model2.add(BatchNormalization())
model2.add(Dense(20, activation='softmax'))
# Compile model
model2.compile(loss='binary_crossentropy', optimizer='adam',metrics=['accuracy'])

# Model 3: 
model3 = Sequential()
model3.add(Dense(246, kernel_initializer="uniform", activation='sigmoid'))
model3.add(Dense(64, activation='sigmoid'))
model3.add(Dense(32, activation='sigmoid'))
model3.add(Dense(20, activation='sigmoid'))
model3.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])

# After comparison we use the first model: 
model.fit(x=X_train, 
          y=y_train, 
          epochs=20,
          batch_size=128,
          validation_data=(X_test, y_test), 
          )

test = pd.read_csv('.../roc.csv')
testset = test.loc[:,'B02001_001E':'political_affiliation_p111975_n']

model.fit(X_stan, y)
predictions_nn = model.predict(testset)

output_nn = pd.DataFrame({
                       '14608': predictions_nn[0],
                       '14604': predictions_nn[1],
                       '14614': predictions_nn[2]},index=y.columns)

output_nn['mean']=np.mean(output_nn,axis=1)

output_nn.sort_values(by='mean',ascending=False,inplace=True)

output_nn

label2 = raw.loc[:,'GEOWALG']
y2 =  label2*1
x_train, x_test, y2_train, y2_test = train_test_split(X_stan, y2, test_size=0.2, random_state=101)

def plotROC(model):
    """
    1. Plot ROC AUC
    2. Return the best threshold
    """
    model.fit(x_train, y2_train)
    predictions = model.predict(x_test)
    probs = model.predict_proba(x_test)
    preds = probs[:,1]
    fpr, tpr, threshold = roc_curve(y2_test, preds)
    roc_auc = auc(fpr, tpr)
    
    # Plot ROC AUC
    plt.title('Receiver Operating Characteristic')
    plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
    plt.legend(loc = 'lower right')
    plt.plot([0, 1], [0, 1],'r--')
    plt.xlim([0, 1])
    plt.ylim([0, 1])
    plt.ylabel('True Positive Rate')
    plt.xlabel('False Positive Rate')
    plt.show()
    
    # report
    acc = accuracy_score(y2_test,predictions)
    print(classification_report(y2_test, predictions))

    # Find optimal threshold
    rocDf = pd.DataFrame({'fpr': fpr, 'tpr':tpr, 'threshold':threshold})
    rocDf['tpr - fpr'] = rocDf.tpr - rocDf.fpr
    optimalThreshold = rocDf.threshold[rocDf['tpr - fpr'].idxmax()]
    
    return acc

rfc2 = RandomForestClassifier(n_estimators=500)
print("accuracy_rfc: ", plotROC(rfc2))

from lightgbm import LGBMClassifier

lgb = LGBMClassifier(objective='binary',
                      learning_rate=0.049,
                      n_estimators=1500,
                      num_leaves=8,
                      min_data_in_leaf=4,
                      max_depth=3,
                      max_bin=41,
                      bagging_fraction=0.845,
                      bagging_freq=5,
                      feature_fraction=0.24,
                      feature_fraction_seed=9,
                      bagging_seed=9,
                      min_sum_hessian_in_leaf=11)
print("accuracy_lgb: ", plotROC(lgb))

lgb.fit(X_stan, y2)
predictions = lgb.predict_proba(testset)[:, 1]

predictions

output = pd.DataFrame({'zipcode': test.zipcode,
                       'probability': predictions})

output